function Rd = Reflectance_Diffuse(r,a,lambda)
% r: roughness, unit (m), 0 < r < 2e-7
% m: slope, dimensionless, 0 < m < ?
% a: auto-correlation length, unit (
% lambda: wavelength, unit (m)

if nargin == 2
    lambda = 550e-9; %550 nm
end
% a = sqrt(2)*r/m;
k = 2*pi^4*(a/lambda)^2*(r/lambda)^2;
Rd = quad(@integrant,0,pi/2);
if Rd>1
    Rd = 1;
end
    function y = integrant(theta)
        y = k*(cos(theta)+1).^4.*sin(theta).*exp(-((pi*a*sin(theta)/lambda).^2));
    end

end